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ABSTRACT 

A simplified non-linear numerical model for the development of incompressible magnetohydrodynamics 
(MHD) in the presence of a strong magnetic field B|| and stratification, nicknamed Shell-Atm, is presented. 
In planes orthogonal to the mean field, the non-linear incompressible dynamics is replaced by 2D shell-models 
for the complex variables u and b, allowing one to reach large Reynolds numbers while at the same time 
carrying out sufficiently long time integrations to obtain a good statistics at moderate computational cost. The 
shell-models of different planes are coupled by Alfven waves propagating along B || . The model may be applied 
to open or closed magnetic field configurations where the axial field dominates and the plasma pressure is low; 
here we apply it to the specific case of a magnetic loop of the solar corona heated via turbulence driven by 
photospheric motions, and we use statistics for its analysis. The Alfven waves interact non-linearly and form 
turbulent spectra in the directions perpendicular and, via propagation, also parallel to the mean field. A heating 
function is obtained, and is shown to be intermittent; the average heating is consistent with values required 
for sustaining a hot corona, and is proportional to the aspect ratio of the loop to the power -L5; characteristic 
properties of heating events are distributed as power-laws. Cross-correlations show a delay of dissipation 
compared to energy content. 

Subject headings: Sun : corona, flares - MHD - turbulence 



L INTRODUCTION 

MHD turbulence in the presence of a mean magnetic 
field, with or without density and gravitational gradients, 
plays a role in man y environm e nts, ranging from stellar 
coronae and w inds (Klein et al. 1991) to the interstellar 
medium (|Desai et al. 1994) and accretion disks. In such re- 
gions, energy may be transferred, accumul ated and dissi- 

■ pated i n a w ay which is inher ently anisotropic ('Shebal in et al.l 

■ 1 19831; lOug hton et a l] 119941; K innev & Mc Williams 119981; 
; iMuUer et al.li2003c lOughton et al.ii2004<) . 

In particular, in solar coronal physics, where one of 
the main problems is to understand how the corona can 

' be sustained at more than a million of degrees, it is be- 
lieved that the necessary heating could be produced at small 
scales generated by a non-linear cascade along a turbulent 
spect rum (Hevvaerts & Priest 1992; Gomez & Ferro Fontan 

' 1 19921) . Furthermore, as flux tubes (e.g. in the form of coro- 
nal loops or coronal funnels) are omnipresent, the anisotropy 
coming from the dominant magnetic field may be a central 
feature of the processes governing energy dissipation, like 
the non-linear collisions of counter-propagating Alfven wave 
packets. It can thus be expected that solving the coronal heat- 
ing problem, i.e. understanding how the temperature of the 
corona can be sustained, may require to understand the details 
of the turbulent dynamics of MHD in these environments. 

One way to study and the dynamics of such system is to 
perform direct numerical simulations (DNS). In the case of 
anisotropic MHD, DNS have provi ded insight on subjects 
like the anisotropy of the sp ectra (e.g. iKinnev & McWilliamsl 
119981: [Oughton et al."2004), the parametric decay of Alfven 
waves (e.g. Del Zanna et al. 2001), and Alfven waves fila- 
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mentation (e.g. iPassot & Sulemll2003h . MHD simulations 
are also used to study the topology of magnetic fiel d lines 
and m agnetic reconnection in the corona (e.g. Aulanie Fet al.l 
l2005h . But the Reynolds numbers attained in all the DNS up 
to now are below 10^, while they are believed to be 10'^ to 
10'"^ in the corona. DNS are very far from being able to repre- 
sent all the scales of turbulence in the corona, there is a huge 
gap to fill. Furthermore, as statistics are of great help in the 
study of turbulence, attempts have been made to analyze sta- 
tistically energy dissipations produced by DNS. Distri b utions 
of e vents are for instance p resented in iDmitruk et al.l (Il998h 
andlGeo rgoulis et al.l (Il998i) from 2D DNS of reduced MHD. 
But it is still difficult to get significant statistics from 3D 
DNS, and it is even more difficult when trying to go to higher 
Reynolds numbers because then the grid resolution must be 
higher and the computations of the model are too slow. For all 
these reasons, there is a need for simplified numerical models 
of MHD, which would run sufficiently fast to get statistics of 
turbulence at high Reynolds numbers, while keeping the most 
relevant features of MHD turbulence. 

Several approaches have been used to build such simplified 
numerical models of MHD. For example, the Self-Organized 
Criticality behavior (SOC) of MHD systems can be modeled 
by cellular automata (CA), where the interactions of indi- 
vidual cells translate into a global statistic behavior of the 
whole system, followin g the first models of iLu & HamiltonI 
(119911) ; iLu et al.1 (Il993h . However, the need for physical re- 
alism is not entirely addressed by the cellular automata, de- 
spite efforts tryi ng to include the constraints issued from the 
MHD equations dVlahos et al.| [T99l llshker et aL.2000„ .20011 
iBuchlin et al.ll2003h . 

Another approach is to simplify the non-linear interactions 
by reducing the number of modes which are allowed to in- 
teract non-linearly. In the coron al loop context, a shell- 
model approach has been used by Nigro et aTl (|2004 120051) . 
We have developed a similar numerical model independently, 
starting from the reduced MHD equations but allowing for 
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the stratification of the plasma. This numerical code, nick- 
named Shell-Atm, allows one to reach (kinetic and mag- 
netic) Reynolds numbers unachieved before. In this paper, 
we will focus on the problem of a coronal loop where energy 
is forced into the system by footpoint motions, describing in 
detail the dynamics of the heating events, turbulence spectra, 
statistics and scaling laws. 

2. DESCRIPTION OF THE Shell-Atm MODEL 

We start from an approximation to incompressible MHD 
known as reduced MHD (RMHD: iKadomtsev & Pogutsi 
Il974t IStraussI [19761) . which is valid when the plasma (3 pa- 
rameter (kinetic over magnetic pressure) is low, the domain 
has a large aspect ratio (a = l/L<^ 1) and the poloidal field 
is small compared to a strong axial external B\\ magnetic field 
(B±/B^\ < a). In this approximation the largest extension L of 
the domain defines the parallel direction, or z-axis; the veloc- 
ity field is only composed of fluctuations u± orthogonal to the 
z-axis; the magnetic field can be decomposed into B|| +B±, 
where B|| = B^\e. is the average magnetic field, parallel to 

the z-axis, and B± is a perpendicular fluctuation. Through- 
out, we will normalize the magnetic fields by y/^opo, con- 
sidering for the moment a medium with uniform density po 

(/7|| = B||/^/MO/Oo and b± = Bj_/^/iopo)- The equations of 
RMHD become: 



du± 



+M I • Vm I =-V 



1- -b , 

PO 2 ^ 



db± 



= b\ ■ \7u\ 



+b 



+b 



dbj 



= 



du± 
V • ZJ_L = 



(1) 

(2) 
(3) 



As one can see in these equations, the non-linear dynamics 
occurs only in the planes perpendicular to the mean field B\\ 
while Alfven waves propagate along the mea n field. Direct 
simulations of these equations in one plane ( Dmitruk et al 



simulations or these equatio ns in one plane ( UmitruK et al 
19981: iGeorgouhs et al.llI998h or in a 3D box fomitruk et al 



20031) have been carried out but the Reynolds numbers ob- 



tained with such simulations are much too low to get a real- 
istic inertial range of turbulence and long-term statistics. It 
is therefore our interest to simplify this model further by re- 
ducing the dynamics in the planes. This can be done by using 
shell-models, as described below. 

The plasma of the solar corona and solar wind is stratified, 
so that one must allow for gradients of the mass density p 
even while considering incompressible couplings. Stratifi- 
cation couples incompressible Alfven waves by introducing 
variations in Alfven speed and therefore reflection (as well as 
amplification/depression of amplitudes due to the conserva- 
tion of energy flux). Such terms may be written more clearly 
in terms of the the Elsasser variables = u±±b± (with 
b± = B± I y/pop), in which case the effect of stratification 
on th e linear propagation of modes may be written as dVellil 
[19931) : 



The first two terms describe the wave propagation, the third 
term the reflection of the waves by the perpendicular gradi- 
ent of the Alfven speed (which vanishes for a non-diverging 
flux tube), the fourth term describes the growth or decrease 
in the normalized wave amplitude due to variations in Alfven 
speed - assuring conservation of wave energy flux - as well as 
the isotropic part of the reflection. We will incorporate these 
terms in the general framework of Eqs. ([TJll), but first we dis- 
cuss how the non-linear couplings are modeled in the shell 
approximation. 

2. 1 . Classical MHD shell-models 

In shell-models of incompressible MHD turbulence 
(Gloaguen et al. 1985; Biskamp 1994; Giuliani & Carboni 
11998; Boffetta et al. 1999; iGiuliani et al. 2002), one starts by 
taking the Fourier transform of the MHD equations and di- 
viding wavevector space into concentric shells S„ = {k\ \\k\\ e 
[kn,k„+i]} with k„ = k()X", n = 0,...,nj^- 1, and usually 
A = 2. Also, a single complex scalar value u„ is cho- 
sen to represent the original longitudinal velocity increments 

(^{x+ £) — u(xij ■ £/i on scales £ for 2tt /£€ S„. The same ap- 
proximation is made for the magnetic field: a scalar value /?„ 
represents the magnetic field increments on the same scales £. 
In this way the non-linear interactions, originally a vector con- 
volution in the 3D vector space are reduced to to a ID summa- 
tion in terms of the shell index n. This one-dimensional model 
is the magnetohyd rodynamic analog of the GOY (Gledzer- 
Ohkitani-Yamada: lYamada & Ohkitani|[l 988) shell-model of 
fluid turbulence. 

One obtains the fo llowing equation, given in 
iGiuliani & Carbond ( [19981) : 

dZ.± 



df 



kl(iy^Z^ + iy-Z^) + ik„T,^*+f-- 



(5) 



where = u„ ± b„ are the Elsasser-like variables, ly^ = 
{v±r])/2 are combinations of kinematic viscosity and resis- 
tivity, are external driving forces, and is the non-linear 
term, obtained by assuming (1) that the non-linear interac- 
tions occur in triads of neighboring modes and (2) the conser- 
vations of the total pseudo-energies = |Z,^ p (and thus 
the energy E = E^ + E~ and the cross-helicity he = E^ - E~) 
and a third invariant H'^ = ^^sign((5- l)"A;"|v„p which de- 
pends on the dimensional ity of the MHD system to model 
(iGiuUani & Carbondll998l) . 



:0 (4) 



2.2. Specificities of the SHELL- ATM model 

The "classical" GOY-like shell-model that we have just pre- 
sented corresponds to MHD, where the average magnetic field 
Z?|| has not yet been separated out; in the SHELL- ATM model 

we present now in this paper, the average magnetic field B\\ 
is separated out, by starting from the RMHD equations ([HO. 
The new model we obtain corresponds basically to a pile of 
planes coupled by Alfven waves and containing each a "clas- 
sical" shell-model for 2D MHD (Fig. [T] top). This is simi- 
lar to t he loop model developed independently bv lNigro et alj 
(|2004|) . but the stratification of the atmosphere that we intro- 
duce allows to use this model in a large variety of cases of 
coronal loops or other structures (although we do not use this 
specific feature in the runs presented in this paper). In the 
Shell-Atm model: 

• The profile of the Alfven speed b\\{z) along the mean 
field (i.e. the atmospheric structuring of the plasma) is 
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given via a density stratification po{z), an areal expan- 
sion factor of the flux tube A{z) and magnetic flux con- 
servation. The latter two effects imply that the width of 
the loop and corresponding wavenumber ko must also, 
in general, depend on z. 

• The Elsasser variables Z,^ now depend on the position z 
of the plane along the main axis of the simulation box, 
and the left hand-side of Eq. (|5]l is replaced by the term 
(d, ± b\\d:)Z^ ± |Z±a(lnp) ± \Zfd-b\y correspond- 
ing to the linear Alfven wave propagation in a strati- 
fied static atmosphere (Eq.|4]i. As a result, the external 
driving forces are not needed anymore, as energy 
can simply be input as an incoming energy flux at the 
boundaries. 

• The non-linear interactions occur inside each plane, in 
2D. In this case the third invariant of MHD is anastro- 
phy, i.e. the total square module of the magnetic po- 
tential {H^ with a = 2). The coefficients of the non- 
linear terms of flie sh ell-model are then those of 
iGiuliani & Carbon^ (1 19981) with parameters a = 2 and 
5 > 1 (i.e. (5 = 5/4 and (5„ = -1 /3). 

To summarize, the fields of the Shell-Atm model we 
introduce and use in this paper are the complex variables 
Z^{z,t), which are the Elsasser-like fields u„(z,t)±b„{z,t). n 
is is the index of the shell, corresponding to the perpendicular 
wavenumber k„{z) = A:o(z)A", with A = 2; it can be any inte- 
ger (positive or negative), but for numerical computations it is 
convenient to assume that Zf{z,t) is outside some domain 
[[0,«j^- 1]]"^. z is the position on the axis supporting in a 
domain [0,L] discretized over planes. The equations of the 
model are: 

{d, ±b\\d,)Zt ± \ztdsnp) ± \z'^d,b\\ = 

-kl(iy^zt + iy-Z^) + ik„T,^* (6) 

with given by IGiuliani & CarbondlT998l (with a = 2, 6 = 

5/4 and5;„ = -l/3). 

2.3. Quantities derived from the fields of the model 

As Z^(z,t) represents the Elsasser field at perpendicu- 
lar wavenumbers included in the shell S„ and at position z, 
\Z^(z,t)\^ /A is the energy per unit mass at position z at time 
f in the modes included in shell S„. If we assume that the 
modelled loop is a cylinder of diameter 2Ti/kQ and that, after 
discretization in the z direction, the separation between planes 
(i.e. the thickness of each plane) is 5z, then the cross-section 
of the loop by a plane is A = 7r^/A:Q and the volume associated 
to each plane is V = A5z\ with a mass density po, the mass as- 
sociated to each plane is m = poAfe and the energy contained 
in the field Z^{z,t) is: 

£±(z,f)= im|Z±(z,f)P= ^Po^&|Z±(z,f)P (7) 

E„ as a function of n (for any field, position and time) 
will hereafter be referred to as the "shell energy spectrum". 
To get a ID perpendicular spectrum (as those given by tur- 
bulence theories), we need in addition to take into account 

The energy flux from tiie domain JO, ~ ll outwards is then zero, as 
can easily be seen using the equation of the spectral energy flux O for n = 
and n = nx ■ 



the geometry of the shell S„ in Fourier space: for a shell- 
model representing 2D MHD, S„ has an area S{S„) = nk^iX^- 
1). It follows that the 2D energy spectral density in the 
shell is E„/S{S„) and that the ID energy spectral density is 
2£'„/(A:„(A^- 1)). Note that for this reason there is adifference 
of 1 between the slope of a power -law ID perpendicu lar spec- 
trum (e.g. -5/3 for a spectrum as lKolmogorovl[T94Tl) and the 
slope of its "shell energy spectrum" counterpart (e.g. -2/3). 

2.4. Scales of quantities of the model and time scales 

The equations are rendered non-dimensional introducing 
characteristic units of time, length and density. For the coro- 
nal situation we choose 10^ m for the unit of length, 1 s for 
the unit of time and 10^ kg for the unit of mass. Then the 
units of the other quantities derive naturally from these basic 
units and are 10 Mm • s~' for velocity, 10"'^ kg • m"-' for mass 
density, 10''*m^ • s~' for diffusivities, 10^^ J for energies, and 
10^^ W for powers. 

The characteristic time scales for each of the terms of 
Eq. (|6]l are deduced from their orders of magnitude: 

d,Z^b\\^Z^Dk\z^k^Z^ (8) 

where Z is an order of magnitude of the fields Zf- of the model 
at wavenumbers k\\ (parallel to S||) and k± = koX", and where 
D represents either ly or rj. We obtain: 

• the Alfven time ta = 27r/(Z7||A:||); 

• the characteristic time of dissipation = (Pk^)~^; 

• the characteristic time of non-linear interactions tnl = 
{k±Z{k±))~^ in the planes. 

• the wave reflection time scale f« = 2/dJb\\. 

The maximum Alfven time r^^max is obtained for k\\ =2'k/L 
and corresponds to the time needed by the wave to cross 
the simulation box. By taking on the other hand the min- 
imum of all the characteristic times in the box (using k\\ = 
27r/& to get rA,mi„, k± = k„^-i to get r^^min, and TNL.min = 
((A:^Z(fc^)),nax)~'), we can estimate the time step needed by 
the numerical schemes according to Courant-Friedrichs-Lewy 
(CFL) conditions. 

Other time scales also appear in the Shell- ATM model: 

• the correlation time t* of the forcing, which depends 
on the precise form of the forcing (see Sect. l2.5l for the 
case of a coronal loop); 

• the turbulent cascade time scale Tcascade = J2n ''nl(^± = 
k„) where the sum is taken on the modes n of the inertial 
range of the spectrum (see Sect. l3.2l l. 

2.5. Case of a coronal loop 

Geometry — For the case of coronal loops forced via photo- 
spheric motions, we consider the loops to be "straightened 
out" as seen in Fig. [1] (bottom). T his is similar to the c ellu- 
lar automaton model presented in iBuchlin et al.l (120031) . but 
here the non-linear interactions between modes of MHD tur- 
bulence are represented through shell-models instead of cellu- 
lar automata. Furthermore, for simplicity, we choose uniform 
density po, Alfven speed b\\ and loop width (27r/A:o); more 
realistic cases wiU be studied in future works. 
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Fig. 1 . — Top: layout of the Shell-Atm model in the general case; shell- 
models in planes orthogonal to iSy are piled up along B|| . Bottom: in the case 
of a coronal loop, the loop is unbent into the cylindric simulation box. 

Forcing — With this geometry, the boundary planes of the 
model represent the footpoints of the loop, which are an- 
chored in the photosphere. We choose to impose a time- 
dependent vortex-like velocity field on modes at the larger 
scales, corresponding to photospheric convective motions at 
the scale of the supergranulation. Since the velocity is im- 
posed, waves travelling along the loop are partially reflected 
at the photosphere. 

The imposed velocity field u„(z, t) on each mode n of both 
boundary planes z = and z = L has the form; 

u,„{t) = ufn (e^'^'^-' sin2(7rf /f *) + e^'''^- sin^CTrf /f * + tt /2)) 

(9) 

where u / „ is the amplitude of the forcing (non-zero only for 
some n corresponding to scales 27r/fc„ of the supergranula- 
tion), and A. „ and B, „ are each independent random complex 
coefficients of module 1 and whose complex argument is uni- 
formly distributed over [0,27r]. These coefficients are kept 
constant during time intervals of duration f *, and they are ran- 
domly changed when t = [t*] and t = t* /2 [f *] respectively, 
i.e. when the corresponding si n^ term is zero. Th is is another 
difference with the model of Nig ro et aP (l2004l) . who force 
using a stochastic velocity function on one boundary plane 
only. The auto-correlation time of the forcing field is then of 
the order off*, which is chosen to be much longer than all the 
other time scales of the model (Sect. 12.4b 

In practice, this boundary condition is realized by impos- 
ing an incoming Alfven wave Z^(0,f) = -Z~{0,t) + 2u() „{t) at 
the boundary z = and an incoming Alfven wave Z~{L, t) = 
-Z^(L, t)+2ui,n(t) at the boundary z = L. The resulting power 
entering the loop is: 

^/ = ^ E (|z:(0, 01' - |z,7(0, 01') 

n 

+i ^ po(L)^|| (i)A(L) {\Z;,{L, 01' - |z:(i, 01') (10) 

n 

Note that this power is not imposed, but it depends on the 
fields already contained in the simulation box. 

3. NUMERICAL SIMULATIONS AND ANALYSIS OF THE RESULTS 




MOO JOZO iW WfiC 3W9 310!) 



imo iv/si JO 40 jaeo 5oat) noo 

Fig. 2. — Energy balance in the model. Top: energies and integrated 
dissipation powers (top thin line: integrated power of forcing; bottom thin 
line: integrated dissipation power; thick line: energy, and sum of integrated 
contributions of powers). The small deviation (only 1 % over the time span 
of this plot) of the energy compared to the sum of integrated powers shows 
that the numerical dissipation is low. Bottom: power time series (from top 
top bottom: forcing power, numerical dissipation power, dissipation power). 
Quantities with negative contributions to the energy balance are shown as 
negative, (color version online) 

3.1. Energy 
3.1.1. Energy balance 

Alfven wave propagation as well as the non-linear terms in 
our shell-model conserve the energy, so that changes in total 
energy in the loop arise only from flux through the photo- 
spheric boundaries (i.e. the forcing) and from the dissipation. 
This energy balance is well verified in practice, within 1 % in 
general as can be seen in Fig.|2] as long as the numerical dis- 
sipation due to the numerical scheme for wave propagation 
is not too high; the condition for this is that the perpendicular 
dissipation scales are not too small compared to the separation 
between planes in the z direction. 

3. 1.2. Ratio of magnetic over kinetic energy 

The ratio of the magnetic to kinetic energy in the station- 
ary state may be estimated as follows. First, a simple linear 
estimate of the velocity field leads to: 

M± = u^fi{x,y)co?,{ujpht)z/L, (11) 

while the magnetic field is given by: 

b^ = u^fi{x,y)b\\t/L. (12) 

The relative importance of the higher freque ncy modes to 
this low-frequency energy flux was discussed bv' Milano et alj 
(|I997): given that the Alfven wave travel times along loop 
is of the order of seconds, while most of the power in pho- 
tospheric motions is in the minutes to hour ranges, it is the 
low-frequency resonance which plays a more important role. 
Energy injection from the photosphere into the corona there- 
fore grows as f^, and is stored in the transverse coronal mag- 
netic field, while the velocity field is bounded by its pho- 
tospheric value. The linear solution will eventually break 
down because the magnetic field determined in Eq. ( fT2] i is 
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= 0.1, Bo=1, Uf=10- 



Fig. 3. — Ratio of average magnetic over kinetic energy as a functio n of 
aspect ratio, plotted witli tlieoretical scaling (plain line) derived from Eq. (13) 
and power-law fit, (dashed line, with slope -0-74 it 0.L4). 

not in general rorce-rree and therefore will cause the coro- 
nal field to evolve dynamically. The ratio of magnetic to ki- 
netic energies at this point may be estimated dimensionally 
by asking for the change in coronal velocity field determined 
by nonlinear interactions in Eq. ([T]i to be of the same order 
of magnitude as the field given by Eq. (fTTT i. Denoting the 
rms photospheric speed by «/, after a time Af the nonlin- 
ear term has the dimensional value u^b^^At^/iL^, growing 
quadratically with time. It will cause a change in the coro- 
nal loop velocity field of the same order as m/ over the time 



At, when u^blAt^/iL^ ^ Uf/ At. One then recovers the time 



2/3 1/3 

Af as Af ~ T^ mnx'^e wherc r^^max is the loop crossing time 
while Te = £/ufis the non-linear time calculated on the photo- 
spheric velocity. The ratio of magnetic to kinetic energies in 
the corona at this stage is then 



^11 



At 



= 3 



TA.i 



2/3 



(13) 



which we associate with the saturation level of magnetic to 
kinetic energies. This ratio should be a function of the as- 
pect ratio a = b^^TA.m-dx/i of the loop. To check whether the 
Shell-Atm model follows this dependence, we perform a 
series of simulations with parameters b\\ = 1, ko = 2Qtt (i.e. a 
width i = 0.1), = ?7 = lO"'', and m/.„ = 10"^ for n £ P,4]]. 
The number of shells is nj^ = 16 and the number of planes 
n, is taken in the set {200,500, 1000,2000,5000} with a sep- 
aration 10"^ between planes in all cases, leading to lengths 
L £ {0.2,0.5, 1,2,5} and aspect ratios a e {2,5, 10,20,50}. 

The ratio of magnetic energy to kinetic energy in the sta- 
tionary state is plotted as a function of aspect ratio in Fig. [3] 
together with what is expected from Eq. ( fT3] l. The numerical 
results we obtain support roughly this scaling, although the 
experimental ratios are smaller than the theoretical ratios by 
a factor 1.36 and the values for an aspect ratio a = 2 deviate 
from the scaling obtained for other aspect ratios. The slight 
departure from the proposed scaling, and the fact that the sat- 
uration level of magnetic to kinetic energy is lower than pre- 
dicted, could be due to the "leakage" of energy to the higher 
frequency resonances of the loop (as shown in Fig.[T3li. 

3.2. Spectra 

3.2.1. Formation of the spectra and spectral energy flux 

The energy flux in each plane from the shells k < k„to the 
shells k> k„ is the derivative of the energy contained in the 
shells k > k„ due to the non-linear terms r„ , namely: 



n„ = 



.s=±l 



+(2-S-S,„)Z'„_2Z„liZ', 

+\{(S+s,„)zi,z:,zz, 

+{2-5-5,„)ZUZ7K,,] 



(14) 
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With ^„ = we re cover the hydrodynamic s pectral energy flux 
given in Eq. (2) of Fr ick & SokolofjjlQQS !). Furthermore, this 
energy flux is consistent with the general idea that the energy 
flows "downhill" in the ID perpendicular energy spectrum. 

When starting the simulation from a very low amplitude 
field (Fig.m, the magnetic energy and then the kinetic energy 
at the scales of the forcing grow, and when the fields are suf- 
ficiently large, non-linear effects become visible as energy is 
transferred to modes beyond those initially forced. In partic- 
ular, there is a flux to higher k±^ mode numbers (11,, > 0, di- 
rect cascade), which continues to highest wavenumbers where 
dissipation occurs, as well as a flux to smaller wavenum- 
bers (n„ < 0, inverse cascade) which energizes modes at the 
largest scales and saturates at a level comparable to the forced 
modes. 

As a result of the energy cascade, an inertial range ap- 
pears between the forcing and the dissipation scales, in the 
same way than in the original shell-models. The energy 
flux n„ across shells is uniform on average over the whole 
inertial range. The Reynolds number in the case shown 
here can be evaluated to 10^, which is much higher than 
any Reynolds number of direct numerical simulations. Even 
higher Reynolds numbers can be attained using more shells 
and planes, at the cost of the ability to do long-term statistics. 

3.2.2. Fhictuations of the spectra 

Once in the stationary phase, the spectra continue to fluctu- 
ate, with characteristic time scales linked to the "local" time 
scales, i.e. the time scales, as described in Sect. 12. 4[ consid- 
ered as depending on the mode k„ (Fig.|5]). The most relevant 
time scale seems to be the local non-linear time scale T^h{k), 
defined in |2.4| as this is the time scale on which the energy in 
a given mode n can change under the action of the non-linear 
terms of Eq. (|6]l. More precisely, no dynamics occur at time 
scales below the local non-linear time scales, as can be seen 
in Fig. |6l the modes with low k±^ have thus only long fluc- 
tuation times, excited by the long time scales of the forcing, 
while the modes with high k±^ fluctuate fast, but still with long 
characteristic times due to the flow of energy coming from the 
modes at low k±^ (these long-term fluctuations are common to 
the whole spectrum). 

3.2.3. Evolution of the spectra during an event. 

To understand what happens during episodes of high energy 
dissipation, we have analyzed the spectra of the fields before, 
during and after such an episode. The differences of spectra 
in respect to an average spectrum (Fig. |7]i show that before 
the event (the maximum dissipation corresponds to the red 
spectrum differences), the energy accumulates over the whole 
spectrum. The total energy is then high, the non-linear times 
are short, and the energy flows down rapidly to the smallest 
scales according to Eq. ( fT4b : it enhances the spectra at the 
largest wavenumbers by several orders of magnitude, leading 
to a strong enhancement of the dissipation power. As energy 
is released, this process leads then to a decrease of the spec- 
trum, first in the dissipative range (high wavenumbers), then 
in the whole spectrum. The dissipation power is then low 
again, and as the non-linear time scales in the inertial range 
are longer, the energy injected at the largest scales cannot flow 
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Fig. 4. — Kinetic (left) and magnetic (riglit) perpendicular spectra of energy in the sliells of tlie model, averaged over the length of the loop. On each plot, 
40 spectra are shown, separated by 10"^ units of time, starting shortly after the beginning of the simulation (lowest curves). The forcing is performed on modes 
con'esponding to logjQA:,, = 2.4 to 3.0. 
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Fig. 5. — Time scales as a function of k^: Alfven time T^r 
linear time t-nl, dissipation time t^, forcing correlation time t* 
time T/i.max (from top to bottom in the caption inset). 
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Fig. 6. — Time series of the total energy contained in the shells n £ 
{2,8,14,20} of the model, and Morlet wavelet time-scale planes of each 
of these time series. The axes of the time-scale planes are time (horizon- 
tally) and time scale (vertically, logarithmic). The theoretical non-linear time 
scale Tm^(k„ ) as a function of time has been superimposed on each time-scale 
plane, (color version online) 

to the smallest scales as fast as before: the energy does not 
reach easily the dissipative scales and the dissipation power 
remai ns low, u ntil the next such episode. 

Nig ro et all (|2005) underline that the leading term of the 
energy flux across scales (Eq. [T4l i is proportionnal to k„blu„ 
(with the notations of the shell-model variables) and they also 
observe short-term variations of the kinetic energy spectrum 




Fig. 7. — Top left: differences (in logjQ-space) between kinetic perpendic- 
ular spectra of energy in the shells of the model and their average (in logjQ- 
space). The spectra are averaged over the length of the loop, and are plotted 
at times surrounding an event of dissipation power (coiTesponding to the red 
spectra): 10 spectra are shown before the event, and 10 spectra are shown 
after; the difference spectra are each separated by 10 units of time, and are 
stacked from top to bottom, with a shift of 1 unit of the y-axis between each 
of them. Top right: same plot for the magnetic energy spectra. Bottom: time 
series of energy and dissipation power, with the colors corresponding to the 
time intervals used to compute the spectra. 

around a dissipation event. These variations appear to control 
the energy flux to the smallest scales, and then the dissipation. 
In addition, we have shown that these variations exist on a 
longer term around an event, and that the magnetic energy 
spectrum also varies on the same time scales. The cross-scale 
energy flux may thus be controlled by both the kinetic and the 
magnetic energy spectra. 

3.2.4. Slopes of the spectra 

The slopes of the power-law ID perpendicular spectra of 
the velocity and magnetic field (Fig. |4]i seem to be roughly 
equal on the inertial range, but as the spectra fluctuate with 
time, there are fluctuations of the slopes. The distribution of 
these slopes obtained at different times is shown on Fig. [HJ 
the median slope is -1 .89 (with a standard deviation 0.10) for 
the velocity spectrum and -1.81 (with a standard deviation 
0.13) for the magnetic spectrum. It appears that, on average, 
the kinetic spectrum is slightly steeper (by 4 %) than the mag- 
netic spectrum. If we look specifically at the times when the 
total dissipation power exceeds its 90th percentile, i.e. during 
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Fig. 8. — Distribution function of the slopes of the ID perpendicular spec- 
trum averaged along for the kinetic (plain line) and magnetic (dashed) 
fields. The median slopes are respectively —1.89 and —1.81, and are plotted 
as vertical lines. The distiibution of the ratios between the slopes for kinetic 
and magnetic perpendicular spectra is shown in the inset, together with its 
median 1.044. 

events of energy dissipation, the spectra are slightly shallower, 
with medians of slopes being -1.83 and -1.77 for the veloc- 
ity and magnetic spectrum respectively. This translates the 
fact that more energy is present at small scales during events 
of energy dissipation (red curves of Fig.]?). 

These ID spect ra are different than the ones found by 
iNigro et al] (l2005l) . corresponding to ID spectra of slopes 
-5/3 for velocity and w -3 for magnetic field; an explana- 
tion could be that their inertial range is smaller, and that their 
fitting range includes scales where forcing occurs. 

3.2.5. Parallel and perpendicular spectra 

In this model non-linear interactions occur only in perpen- 
dicular planes. Development of small scales along the mag- 
netic field is then merely a consequence of the Alfvenic prop- 
agation of differences in the dynamics in different planes. 
One therefore expects parallel and perpendicular spectra to 
be different, with a relationship determined by the so-called 
critical balance condition, namely, that for a given perpen- 
dicular scale, differences in the parallel direction can ap- 
pear only between planes such that the Alfven propagation 
time is longer th an the (perpendicular) non- L inear time at that 
same scale (e g. iGoldreich & Sridhad[T99l ICho et aT]l2002l: 
lOughton et al.ll2004l) . In the case of the present model, as- 
suming a k~°' ID energy spectrum (i.e. a A:"""*"' "shell energy 
spectrum"), the non-linear time scale is tnl(^±) oc A:'""^'/^. 
With a constant and uniform advection velocity b\\, the criti- 
cal balance condition can be expressed by: 



Note that with a Kolmogorov a = 5/3 spectrum, we recover 
the result k\\ cx k^l^ of lGoldreich & Sridhaj (11995!) . 

For a field a„{z) of the model at a given time t {a can be Z^, 
u„ or b„), let d„(k\\) be its Fourier transform along the z-axis. 
We obtain the two-dimensional power spectrum of a (function 
of k±^ = k„ and ) from: 

A(k^,k\\) = -^\a„(k\\)f (16) 

where c is a constant. 

To get a sufficient wavenumber range in the parallel and 
perpendicular directions, we need to do simulations with a 
very large number of planes. This is achieved by starting 
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Fig. 9. — Total spectrum of the Z+ and Z fields, as a function of the 
perpendicular and of the parallel wavenumbers. The plain lines are level 

lines, and the dashed lines are ellipses of axes ij^ and oc k'J^' , for different 
values of k± . 

a simulation with a number of planes o, and then, once 
the energy has reached its final order of magnitude, by stop- 
ping the simulation and resuming it after having interpolated 
the fields in the z direction. We can perform several steps 
of this process if needed. Figure |9] shows a 2D spectrum 
obtained by summing the and Z" spectra and by aver- 
aging them over 10 times separated by lOrA^mux, during a 
run with «; = lOn^.o = 5000 planes (runs with 50000 planes 
were also performed). The level lines in the {k±,k^^) space 
are clearly non-circular, and appear to follow the critical bal- 
ance ellipses (Eq. [Tsl ) at large k± though with an excess of 
energy in the paral lel direction. The aniso tropy angle as de- 
fined by Eq. (5) of iDel Zanna et al] (1200 ll) . computed on the 
range logjgA: £ [1.8,4.4] (where the spectrum is known as a 
function of both k± and ^y) is 67 degrees, which confirms that 
the spectrum is elongated in the perpendicular direction. 

3.3. Dissipation, heating function and statistical properties 

Heating function. — If we look at the energy dissipation power 
per unit length as a function of both time and position along 
the loop, we get the "heating function" (Fig. [TOb). We see 
again (and for the same reasons than before) short-lived events 
of dissipations, and they correspond to short structures along 
the axis of the loop, whose size is of the order of the propaga- 
tion distance of the structure during its lifetime. Some Alfven 
wavepackets are also strong enough to be dissipated only af- 
ter interacting with a lot of counter-propagating wavepackets, 
and thus they live a longer time and leave an oblique trace in 
the heating function during their propagation. 

Furthermore, when we look at the heating function at long 
time scales of hundreds of crossing times TA.max (Fig. [Tob). 
some features appear, which are related to the slow variations 
of the total energy (mainly contained in the slow-variating 
low-^^ modes) under the effect of the slow forcing of time 
scale t* (which is chosen to be a few hundreds of TA.max)- The 
time variations of the dissipation power at these times scales 
(corresponding to a few minutes of physical time) seem to be 
almost the same on all positions along the loop. This is con- 
sistent with the common statement that the loop is heated as a 
whole, even though (1) the elementary events of dissipation, 
as seen in Fig. [TOh . are each small compared to the length 
of the loop (2) thermodynamics, which would further smooth 
out the appearance of the heating function obtained from ob- 
servable variables (because of the fast conduction times) has 
not yet been taken into account. 
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Dissipation power time series. — The integral of the heat- 
ing function along the loop is the time series of the power 
of energy dissipation e(f), shown in the bottom panels of 
Figs.[TOl'a-b). These time series display spikes of high dissi- 
pation power at short time scales during high-activity periods, 
as is usually found in both observations of solar flares and 
simulations of high-Reynolds number MHD turbulence. This 
dissipation time series will be analyzed further later. 

Average profile of dissipation power — On the Other hand, the 
time average of the heating function, i.e. the average power of 
energy dissipation per unit length as a function of the position 
along the loop, shown in Fig. [TOl'c). is almost flat and drops 
only near the loop footpoints. This would suggest that coronal 
heating takes place almost uniformly along loops, although 
not at footpoints; however one must bear in mind that these 
simulations do not yet take into account realistic profiles of 
density and magnetic field. 

Intermittency — The increments (5re(f) = e(f + r)-e(f) of the 
time series e(f) at a given time lag r have a distribution 
whose shape depends on the time lag r: in Fig. (TT] the dis- 
tributions of the ^re(f) normalized by their standard devia- 
tion have wider wings for short time lags than for long time 
lags. Hence the time series is intermittent, which is con- 
firmed by the rise of the flatness (fourth normalized struc- 
ture function) F{t) = {\5re{t)\'^)t I {\5T<^{t)\''-)'^ for smafl time 
lags r (inset of Fig.fTTli. This intermittency is a consequence 
of the intermittency that can be observed in the velocity and 
magnetic fields of the model, and which is also predicted 
by m odels such as [Sh e & Leveque (1994) in hydrodynamics 
and iPolitano & Pouquet(,1995) in MHD. It could be a conse- 
quence of the fluctuations of the spectral flux resulting from 
the long-term glob al fluc tuations of the spectrum which have 
been seen in Sect. 13.2.21 The modes with high k±^ are then 
intermittent, and as they participate the most to energy dissi- 
pation, the time series of energy dissipation power is intermit- 
tent. 

Events — Statistics issued from observations often involve the 
detection of events, or structures, from the observe d fields 
jAschwanden et al. 2000; Parnell & Jupp 2000; BuchU n et al] 
120061) . and the distributions of t heir character i stics. Follow- 
ing the "threshold" definition of iBuchlin et al.l (l2005h . with a 
threshold fixed at the average dissipation power, we get the 
distributions shown in Fig.[T2]for the event total energy con- 
tent, the peak power of energy dissipation, the duration of 
events, and the waiting-time between two successive events. 

The distribution of the peak power in events is narrow, as 
a result of the summation of the heating function over the 
whole loop: the local spikes of energy dissipation are hid- 
den by the average dissipation occurring in the whole loop. 
On the other hand the distributions of integrated dissipation 
power (total energy content of events) and of event durations 
are v ery wide. This is par tially due to the threshold definition 
used (i Buchlin et al.l2005l) . in the case of this time series where 
long time scales are superimposed to the shorter time scales 
of energy dissipation in the dissipation range. Furthermore, 
the waiting-times between successive events have also a wide 
power-law dis tribution. However, as discussed extensively in 
iBuchlin et al.l (l2005i) . this result depends on what definition of 
an event is used to extract events from the time series of the 
power of energy dissipation. 

Compared to the distributions of eve nts obtained from the 
loop shell-model of iNigro et al.l (|2004|) . the main difference 
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Fig. 10. — (Top and middle panels) Heating function, or power of energy 
dissipation per unit length as a function of time t and position z along the loop. 
Two different time intervals are shown in sub-figures (the top and middle 
panels). The bottom panels in each sub-figure represents the integral along 
the loop of the heating function (i.e. the total power of energy dissipation as 
a function of time). (Bottom panel) Time average, over I2OOT4 n,ax following 
(top panel), of the heating function as a function of the position along the 
loop. 

is the much steeper slope (-4.9 instead of -1.8) of the dis- 
tribution of the peak power in events. The reason could be 
the summation effect due to the existence of more but smaller 
dissipation events along the loop, because of the higher reso- 
lution we used in this run, both along the loop {n- = 2000 in- 
stead of n~ = 200, allowed by the parallelization of our code) 
and in the perpendicular direction {n^= 16 instead of 11). 

3.4. Frequencies and time correlations 
3.4.1. Frequencies 
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Fig. 1 1 . — Distributions of tlie increments of the energy dissipation power 
time series, for different time lags. Inset: flatness corresponding to these 
dissipation power increments, (color version onhne) 
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Fig. 12. — From left to right and from top to bottom: distributions of peak 
power, total energy content, duration, and waiting times, for events found in 
the time series of energy dissipation power. The dashed line represents one 
event per histogram bar. 

The time series of the kinetic and magnetic energies reveal 
oscillations corresponding to exchange of energy between the 
velocity and magnetic field. These exchanges occur, for ex- 
ample, thanks to the crossing of Z"*" and Z" wavepackets, 
and should have periods which vary depending on the pre- 
cise number and repartition of wavepackets along the loop. 
A characteristic time scale is of course given by the Alfven 
crossing time r^.max, corresponding to the first resonance fre- 
quency /o = 1 /TA.m-ix- Multiples of the Alfven crossing time 
correspond to higher-frequency resonances of a loop in linear 
theory. While the power spectrum of the time series of total 
energy is a power-law of index -2 over more than 4 decades, 
the spectra of magnetic or kinetic energy display peaks corre- 
sponding to these resonances. The spectrum of the time series 
of kinetic energy (Fig.fTsTl fits to a power-law of index -2.5 at 
very low frequencies. The first resonant frequency, together 
with the higher frequency harmonics, appear as peaks overly- 
ing a different, steeper power-law for the higher frequencies, 
shown in the bottom panel of Fig. [T3] The frequencies of the 
peaks correspond well with the integer multiples nfo of the 
fundamental for « > 5, while at lower frequencies they ap- 
pear shifted. This shift, which is absent in a linear simulation 
(realized with the same parameters but without shell-models, 
i.e. with no non-linear interactions), is probably due to anhar- 
monicity intr oduced b y the n on-Unear effects, as shown by 
iMilano et all ([T997h or lNigrol (l2005h . 

An even better understanding of these oscillations may be 
gleaned from a time-frequency analysis via wavelet trans- 
form shown in Fig. [141 there are oscillations which have 
long but finite lifetimes and different frequencies dominantly 
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Fig. 13. — Plain black lines: spectrum of the time series of kinetic energy. 
The bottom panel is a zoom on the high-frequency range of the spectrum, 
where vertical lines represent the harmonics of the first resonance frequency 
/o = 1- Power-law fits and the slopes obtained are superimposed; the hori- 
zontal range of the lines indicates the range of the fits, and they are shifted 
vertically for clarity. The spectrum (in arbitrary units) of the square ampli- 
tude of the forcing function {\uo „(t)\^ from Eq.|9) is superimposed on the top 
panel (dashed grey line). 

around the fundamental harmonic. These oscillating high- 
frequency wavepackets appear to arise in association with 
dissipation bursts, seen in the dissipation power time-series 
(bottom panel of Fig. [T4b . This intermittent rise in the high- 
frequency component of the velocity field may be involved in 
the enhanced non-linear interactions required to generate the 
bursts in pow er, as required in th e flare-driving mechanism 
highlighted in iNigro et al] (l2005h . On the other hand, their 
persistence may be related to excitation by the time-space lo- 
calization of the burst themselves, a sort of post-microflare 
resonant ringing, which might be observable with future high- 
cadence spectroscopic measurements. 

The comparison between the spectra of the forcing func- 
tion at a boundary and of the resulting energy time series (top 
panel of Fig. [TsT l makes it clear that the spectrum of energy 
is not contaminated by the spectrum of the forcing function, 
as the latter only contains very low frequencies, at or below 
1/f *; this may not be th e case with a stochastic forcing func- 
tion as the one used by Nigro et al.l (1200 4). This underlines 
the role of turbulence in providing the high frequencies which 
can resonate in the loop, viewed as a cavity for Alfven waves, 
even in the absence of an external driver at these frequencies. 

3.4.2. Auto-correlations 

The correlation time of the energy time series is a few 
dozens of the loop crossing times by the Alfven waves 
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Fig. 14. — Morlet wavelet time-scale plane of the kinetic energy time 
series (logaiithmic color scale). Oscillations of long but finite lifetime and 
of different frequencies can be seen at time scales (vertical axis) between 0.2 
and 0.6. Left panel: average wavelet spectrum. Bottom panel: time series of 
total dissipation power, (color version online) 

(Fig. [TST i, consistent with the slow evolution of the en- 
ergy that we have already noted. The correlation time of 
the dissipation power time series is shorter, but still longer 
than the wave crossing time, as an effect of the weakness 
of the intermittent non-linear interactions between counter- 
propagation wavepackets and of the global long-term fluctua- 
tions of the spectrum (including the dissipation range) noted 
in Sect. [3221 

3.4.3. Cross-correlations 

A common goal when studying solar flares and space 
weather is to find precursors of flares, so as to get previsions of 
possible solar-terrestrial events (e.g. lAbramenko et al.ll2002t 
lAbramenko 2005; Georgoulis 2005). With this heating model 
we can use the cross-correlations between time series to in- 
vestigate which time series react first, and what kind of ob- 
servations would be h elpful when predic ting flares. We have 
extended the study of iNigro et alJ (12005 ). who show that in 
some events the kinetic energy begins to grow just before the 
start of a dissipation event, by doing a systematic correlation 
study between all energy and dissipation time series (kinetic, 
magnetic and total). Figure [15] shows the cross-correlations 
between the dissipation e and the energy E time series, as well 
as between the dissipation e and the kinetic energy time 
series^. Both cross-correlation functions show delays of the 
dissipations compared to the energy: the dissipation lags of 
approximately 5 time units compared to the total energy, and 
of approximately 0.5 time units compared to the kinetic en- 
ergy (this last res ult is a confirmation of the result obtained in 
iNigro et ani2005 ). Thus diagnostics methods based on the to- 
tal energy or including the magnetic fields measurements may 
provide more useful results for space weather prediction than 
methods based on the velocity field alone. However, in both 
cases the delays involved are short, of the order of a minute at 
best. 

3.5. Parametric study of dissipation power 

Using the same runs as in Sect l3.1.2l (a set of loop with 
different aspect ratios for a fixed width), we compute the av- 
erage energy dissipation in a stationary state, and we plot it 
versus the aspect ratio (Fig. [T6] i: the scaling of the energy 
dissipation power per unit volume is approximately in a"^/^. 
This scaling ca n be compared to the different models listed in 
iMandrini et al.l (i2000.) and corresponds to heating of MHD in 
2 dimensions. 

^ Other correlations of pairs of time series of kinetic, magnetic or total 
energy or dissipations have not been plotted because they are very similar to 
either of the plotted correlations, as e,, <C tfc ~ e and £„ <C ~ £ . 
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Fig. 15. — Auto-correlation functions and cross-correlation functions of 
time series of energy and energy dissipation power. E is energy, £„ is kinetic 
energy and e is the dissipation power. 

As the slope of this power-law scaling is steeper than 
-1, shorter loops are more efficient in dissipation power per 
unit surface. This can be explained by the fact that Alfven 
wavepackets that reflect on the loop footpoints interact more 
frequently in a short loop than in a long loop; as a matter of 
fact, simulations done when varying the Alfven speed show 
that the average dissipation power also increases when the 
Alfven speed increases. Assuming that the physical units of 
the model (see Sect l2.4l are 10 Mm, 5 s and 10''kg, yield- 
ing ^ = 1 Mm, b\\ = 2 Mm • s"' and po = 10"'^ kg • m"^, we 
get dissipation powers per unit surface between 10^ W • m"^ 
for large aspect ratios and lO^* W-m"^ for small aspect ra- 
tios. These values would b e sufficient to heat the quiet corona 
dWithbroe & Novell 1 1977h . Note however that they also de- 
pend on the physical properties b\\ and po that we have as- 
sumed for the loop. Another series of runs has been per- 
formed to explore the influence of b\\ on the heating, and it 
gives es oc b^''', this translates the fact that wavepackets inter- 
act more frequently when the Alfven speed is higher, leading 
to more dissipation. These both fits, combined with a dimen- 
sional analysis on the variables es (dissipation power per unit 
surface), po (mass density), b\\ (Alfven speed), uj (forcing 
speed) and a (aspect ratio), give: 




for the dissipation power per unit surface in S.l. units 
( W-m"^), as a function of the other variables in S.l. units 
(kg • m"-' and m • s"'). 

4. CONCLUSIONS 

We have presented the Sh ell- ATM model, which is a gen- 
eralization of shell-models dGiuliani & Carbon3ll998l) with 
propagation of Alfven waves along a B\\ field, with the fur- 
ther possibility of introducing a longitudinal stratification of 
physical properties (Z?||, mass density, flux tube expansion fac- 
tor). Although the model is simple and includes only sim- 
plified physical processes, it has a very interesting complex 
non-linear dynamics and it is fast enough to get statistics of 
its fields and of the heating it produces: the simplifications 
we made allow to explore other properties than those acces- 
sible to classical direct numerical simulations (DNS). While 
it is not meant to replace and cannot replace DNS, for exam- 
ple because of the lack of three-dimensional information on 
field-line topologies, it partially fills the huge gap between the 
Reynolds numbers in DNS and in the real corona: although 



Shell-model of RMHD turbulence and coronal heating 



11 



10-= p 




10-^^ , - 

1 1 

Aspect ratio 

Fig. 16. — Average power of energy dissipation per unit volume (model 
dimensionless units) vs. aspect ratio a, for a fixed loop width £ = L/a = OA 
and external field ho = I. The power-law fit (plain line) has a slope —1.52. 
The two dashed lines represent a dissipatation power per unit surface of 10^ 
and lO' W ■ nC^ respectively. 

Reynolds numbers reached in the model, of the order of 10^, 
are still lower than the ones expected in the corona, this al- 
ready represents an outstanding progress compared to DNS. 
Furthermore this allows to explore regimes of MHD turbu- 
lence that are not accessible to DNS, for instance this allows 
for intermittency to appear in turbulence while having a com- 
plete description of the spectra and of the non-linear interac- 
tions in the perpendicular direction. On the other hand, having 
even higher Reynolds numbers, and thus smallest scales, may 
require to take into account non-MHD effects, like kinetic ef- 
fects, that can still be neglected in this model. 

The model has been used in this paper in the case of 
a magnetic loop in the solar corona, in which the physi- 
cal properties of the medium (namely the external longitu- 
dinal magnetic field B|| and the mass density) are assumed 
to be uniform along the loop. In this case, and thanks to 
the above-mentioned characteristics of the model, we could 
show that this model loop displays: a dynamics over a very 
wide range of spatial and temporal scales (4 to 5 orders of 
magnitude); spectra which are formed by a local cross-scale 
energy flux, and which have a wide inertial range in either 
direction, perpendicular or parallel to the external magnetic 
field; a clear anisotropy between the parallel and perpendic- 
ular spectra, which could be compatible with the "critical 
balance" phenomenology; a scaling of the average ratio of 
the magnetic energy over the kinetic energy consistent with 
RMHD; a heating function with multiple spatial and tempo- 
ral scales; a flat longitudinal profile of the average dissipation 



power (although this may be dismissed by further simulations, 
with non-uniform physical properties of the medium along the 
loop); a spiky, and statistically intermittent time series of en- 
ergy dissipation power; power-law distributions of the charac- 
teristics (peak energy, total energy, duration, waiting-times) of 
"events" extracted from the time series of dissipation power; 
finite-lifetime packets of resonant frequencies in the time se- 
ries of energy, but whose frequencies are shifted from the 
harmonics of the linear resonant frequencies because of the 
non-linearities; long-range time correlations in time series; a 
delay of the dissipation time series compared to some energy 
time series; an average dissipation power that scales with the 
loop parameters, and that could be sufficient to sustain the 
high coronal temperatures. As discussed in the text, some 
of these results confirm or complete the results of a similar 
model (Nigro et al. 2004, 2005). 

Further directions for the study of the solar corona using 
this model include taking advantage of the possibility of mod- 
elling non-uniform regions to ( 1 ) allow for density gradients 
in a coronal loop, to seek for the preferred locations of coro- 
nal heating, and (2) study a magnetically open region such 
as a coronal hole. In order to get diagnostics that can be 
compared to observations, this heating model can be coupled 
to the thermodynamics of a loop (including the cooling by 
conduction and radiation) upon which forward-modelling of 
coronal spectral lines may be carried out. We also believe that 
this model can be used in other heliospheric and astrophysical 
applications where MHD applies and where there is a strong 
dominant magnetic field (see Sect. lA.3] in Appendix for code 
availability). 
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APPENDIX 
THE NUMERICAL CODE 
Numerical schemes 

The time advancement of the non-linear terms of the shell-model s is performed by a third-order Runge-Kutta scheme. The 
Alfven-wave propagation is done by the Fromm numerical scheme jFrommll 19681) . Finally, the dissipation terms of the shell- 
models can be computed by an implicit scheme. This allows to relax the CFL condition on and thus to fully resolve the 
dissipation range of the spectrum at no further computational cost. 

Parallelization and parallel efficiency 

The Shell- ATM model is parallelized using MPI, by simply distributing the planes over the processors. Communications are 
mainly needed for the propagation of the Alfven waves between the domains corresponding to the different processors, and for 
the output. The resulting parallelization efficiency is good, and is even close to the ideal parallelization efficiency (up to hundreds 
of processors for n~ = lO'*) thanks in particular to effects due to the cache size (when the number of processors grows, the local 
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data get small enough to fit entirely in the level-2 memory cache of each processor). 

Architecture of the code and availability 

The Shell- ATM code is modular and can be adapted to a large variety of physical systems. Different models for the non- 
linearities and different numerical schemes can be chosen. We believe that the code can be useful for the community and we have 
thus released it under the GNU GPL hcence. The code and its manual can be found at 
|http : //www . arcetri .astro . it/~eric/shell- atm/ coded oc/ 
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